#pragma once

#include "vector.h"
#include <algorithm>
#include <functional>

namespace g4 {

//symmetric matrix
template <typename T>
class basic_smatrix;

typedef basic_smatrix<double> smatrix;
typedef basic_smatrix<float> smatrixf;
typedef basic_smatrix<long double> smatrixl;

template <typename T>
class basic_smatrix {
public:
    T xx, xy, xz, xt, yy, yz, yt, zz, zt, tt;
public:
	typedef T * iterator;
	typedef const T * const_iterator;
public:
	T * begin()
		{ return &xx; }
	const T * begin() const
		{ return &xx; }
	T * end()
		{ return &tt + 1; }
	const T * end() const
		{ return &tt + 1; }
	static int size() 
		{ return 10; }
public:
	basic_smatrix()
		{ }
	basic_smatrix(T _xx, T _xy, T _xz, T _xt, T _yy, T _yz, T _yt, T _zz, T _zt, T _tt) :
		xx(_xx), xy(_xy), xz(_xz), xt(_xt), yy(_yy), yz(_yz), yt(_yt), zz(_zz), zt(_zt), tt(_tt)
		{ }
	template <class U> 
	basic_smatrix(const basic_smatrix<U> & v) :
		xx(T(v.xx)), xy(T(v.xy)), xz(T(v.xz)), xt(T(v.xt)), yy(T(v.yy)), yz(T(v.yz)), yt(T(v.yt)), zz(T(v.zz)), zt(T(v.zt)), tt(T(v.tt))
		{ }
	basic_smatrix & identity(T diagval = 1) 
		{ xx = yy = zz = tt = diagval; xy = xz = xt = yz = yt = zt = 0; return * this; }
	basic_smatrix & zero() 
		{ std::fill(begin(), end(), T(0)); return * this; } 
public:
	basic_vector<T> operator ()(const basic_vector<T> & v) const;
public:
	T det() const;
	operator bool() const ///< true if not degenerate
		{ return det() != 0; };
	basic_smatrix<T> inverse() const;
	/// uses precomputed value of determinant
	basic_smatrix<T> inverse(T det) const;
public:
    basic_smatrix & operator += (const basic_smatrix<T> & v)
		{ std::transform(begin(), end(), v.begin(), begin(), std::plus<T>()); return * this; }
    basic_smatrix & operator += (T a) ///< adds identity matrix multiplied on a
		{ xx += a; yy += a; zz += a; tt += a; return * this; }
    basic_smatrix & operator -= (const basic_smatrix<T> & v)
		{ std::transform(begin(), end(), v.begin(), begin(), std::minus<T>()); return * this; }
    basic_smatrix & operator -= (T a)
		{ xx -= a; yy -= a; zz -= a; tt -= a; return * this; }
    basic_smatrix & operator *= (T a)
		{ std::transform(begin(), end(), begin(), std::bind2nd(std::multiplies<T>(), a)); return * this; }
    basic_smatrix & operator /= (T a)
		{ std::transform(begin(), end(), begin(), std::bind2nd(std::divides<T>(), a)); return * this; }
	bool operator == (const basic_smatrix<T> & v) const
		{ return std::equal(begin(), end(), v.begin()); }
	bool operator != (const basic_smatrix<T> & v) const
		{ return !(*this == v); }
};

template <typename T>
T basic_smatrix<T>::det() const {
	return
		xx*yy*zz*tt 
		- (xx*yy*zt*zt + xx*zz*yt*yt + xx*tt*yz*yz + yy*zz*xt*xt + yy*tt*xz*xz + zz*tt*xy*xy)
		+ xy*xy*zt*zt + xz*xz*yt*yt + xt*xt*yz*yz
		+ 2 * (
			(xx*yt*yz*zt + yy*xz*xt*zt + zz*xy*yt*xt + tt*xy*xz*yz)
			- (xy*yt*xz*zt + xy*xt*yz*zt + xz*yt*xt*yz));
}

template <typename T>
inline basic_vector<T> basic_smatrix<T>::operator ()(const basic_vector<T> & v) const {
	return basic_vector<T>(
		xx * v.x + xy * v.y + xz * v.z + xt * v.t,
		xy * v.x + yy * v.y + yz * v.z + yt * v.t, 
		xz * v.x + yz * v.y + zz * v.z + zt * v.t,
		xt * v.x + yt * v.y + zt * v.z + tt * v.t);
}

template <typename T>
inline basic_smatrix<T> basic_smatrix<T>::inverse() const {
	return inverse(det());
}

template <typename T>
inline basic_smatrix<T> basic_smatrix<T>::inverse(T det) const {
	return basic_smatrix<T>(
		yy*zz*tt - yy*zt*zt - tt*yz*yz - zz*yt*yt + 2*yt*yz*zt, //xx
		-xy*zz*tt + xy*zt*zt + xz*yz*tt - yt*xz*zt - xt*yz*zt + yt*xt*zz, //xy
		xy*yz*tt - xy*yt*zt - yy*xz*tt + yy*xt*zt + xz*yt*yt - yt*xt*yz, //xz
		-xy*yz*zt + xy*yt*zz + yy*xz*zt - yy*xt*zz - xz*yz*yt + xt*yz*yz, //xt
		xx*zz*tt - xx*zt*zt - tt*xz*xz - zz*xt*xt + 2*xz*xt*zt, //yy
		-xx*yz*tt + xx*yt*zt + xz*xy*tt - xt*xy*zt - xz*xt*yt + xt*xt*yz, //yz
		xx*yz*zt - xx*yt*zz - xz*xy*zt + xt*xy*zz + xz*xz*yt - xz*xt*yz, //yt
		xx*yy*tt - xx*yt*yt - tt*xy*xy - yy*xt*xt + 2*xy*xt*yt, //zz
		-xx*yy*zt + xx*yz*yt + xy*xy*zt - xy*xz*yt - xt*xy*yz + xt*xz*yy, //zt
		xx*yy*zz - xx*yz*yz - yy*xz*xz - zz*xy*xy + 2*xz*xy*yz  //tt
	) / det;
}

template <typename T>
inline basic_smatrix<T> operator + (const basic_smatrix<T> & a, const basic_smatrix<T> & b) {
	basic_smatrix<T> t(a);
    return t += b;
}

template <typename T>
inline basic_smatrix<T> operator + (const basic_smatrix<T> & v, T a) {
	basic_smatrix<T> t(v);
    return t += a;
}

template <typename T>
inline basic_smatrix<T> operator - (const basic_smatrix<T> & a, const basic_smatrix<T> & b) {
	basic_smatrix<T> t(a);
    return t -= b;
}

template <typename T>
inline basic_smatrix<T> operator - (const basic_smatrix<T> & v, T a) {
	basic_smatrix<T> t(v);
    return t -= a;
}

template <typename T>
inline basic_smatrix<T> operator * (const basic_smatrix<T> & v, T a) {
	basic_smatrix<T> t(v);
    return t *= a;
}

template <typename T>
inline basic_smatrix<T> operator * (T a, const basic_smatrix<T> & v) {
	basic_smatrix<T> t(v);
    return t *= a;
}

template <typename T>
inline basic_vector<T> operator * (const basic_smatrix<T> & m, const basic_vector<T> & v) {
    return m(v);
}

template <typename T>
inline basic_smatrix<T> operator / (const basic_smatrix<T> & v, T a) {
	basic_smatrix<T> t(v);
    return t /= a;
}

///tensor square (a*a)
template <typename T>
inline basic_smatrix<T> tsquare(const basic_vector<T> & a) {
	return basic_smatrix<T>(
		a.x * a.x, //xx
		a.x * a.y, //xy
		a.x * a.z, //xz
		a.x * a.t, //xt
		a.y * a.y, //yy
		a.y * a.z, //yz
		a.y * a.t, //yt
		a.z * a.z, //zz
		a.z * a.t, //zt
		a.t * a.t  //tt
	);
}

///symmetric sum of tensor products (a*b + b*a)
template <typename T>
inline basic_smatrix<T> sproduct(const basic_vector<T> & a, const basic_vector<T> & b) {
	return basic_smatrix<T>(
		2 * a.x * b.x, //xx
		a.x * b.y + a.y * b.x, //xy
		a.x * b.z + a.z * b.x, //xz
		a.x * b.t + a.t * b.x, //xt
		2 * a.y * b.y, //yy
		a.y * b.z + a.z * b.y, //yz
		a.y * b.t + a.t * b.y, //yt
		2 * a.z * b.z, //zz
		a.z * b.t + a.t * b.z, //zt
		2 * a.t * b.t  //tt
	);
}

} //namespace g4

template <typename T>
class squarer<g4::basic_smatrix<T> > {
	T s;
public:
	explicit squarer(const g4::basic_smatrix<T> & v) : 
		s(v.xx*v.xx + v.yy*v.yy + v.zz*v.zz + v.tt*v.tt + 2*(v.xy*v.xy + v.xz*v.xz + v.xt*v.xt + v.yz*v.yz + v.yt*v.yt + v.zt*v.zt)) 
		{ }
	operator T () const { return s; }
};

template <typename T>
class absoluter<g4::basic_smatrix<T> > {
	T a;
public:
	explicit absoluter(const g4::basic_smatrix<T> & v) : a(::sqrt(sqr(v))) { }
	operator T () const { return a; }
};
